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Abstract. We describe the two-dimensional TreePM method 
in this paper. The 2d TreePM code is an accurate and efficient 
technique to carry out large two-dimensional N-body simulations 

^^ I in cosmology. This hybrid code combines the 2d Barnes and Hut 

^P ' Tree method and the 2d Particle-Mesh method. We describe the 

splitting of force between the PM and the Tree parts. We also 

Tjlj- I estimate error in force for a realistic configuration. Finally, we 

O ■ discuss some tests of the code. 
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1. Introduction 



S^ ' It is believed that large-scale structures in the Universe have formed from 



the gravitational amplification of initial seed density perturbations. Evolu- 
tion of density perturbations at scales smaller than the Hubble radius in an 
expanding background can be studied in the Newtonian limit in the matter- 
dominated regime. Linear theory can be used to study the growth of small 
perturbations in density. But in the absence of analytical methods, numer- 
ical simulations are the only tool available for studying clustering in the 
non-linear regime. The last two decades have seen a rapid development of 
techniques of cosmological simulations as well as computing power and the 
results of these simulations have provided valuable insight into the study of 
structure formation. 

A number of attempts has been made over the past decade to model 
the non-linear evolution of constructs like the two-point correlation func- 
tion using certain non-linear scaling relations [16], [21]. In these relations, 
the evolution of the correlation function can be divided into three distinct 
regimes [22] - the linear regime, the intermediate regime and the non-linear 
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regime. Clearly, a large dynamic range is required in any N-body simulation 
in order to address the issue in all the three regimes under consideration. It 
has been pointed out [2], [20] that by simulating a two-dimensional system 
a much higher dynamic range can be achieved as compared to a complete 
three-dimensional simulation with similar computational resources. 

The simplest N-Body method that has been used for studying clustering 
of large scale structure is the Particle Mesh method. It has two elegant 
features in that it provides periodic boundary conditions by default, and 
the force is softened naturally so as to ensure collisionless evolution of the 
particle distribution. However, softening of force done at grid scale implies 
that the force resolution is very poor. In particular, this limits the dynamic 
range over which we can trust the results of the code [10], [1]. 

A completely different approach to the problem of computing force is 
used in the Tree method. In this approach we consider groups of particles 
at a large distance to be a single entity and compute the force due to the 
group rather than sum over individual particles. There are different ways of 
defining a group, but by far the most popular method is that due to Barnes 
and Hut [6]. 

Several attempts have been made to combine the high resolution of a 
Tree code with the natural inclusion of periodic boundary conditions in a 
PM code [25, 9, 12]. The TreePM method in three dimensions is a hybrid 
N-body method which attempts to combine the same features [3] . The basic 
motivation for these codes is to improve the acceptable dynamic range of 
simulations without a proportionate increase in computational requirements. 
In this paper, we describe the TreePM method for a two-dimensional system. 

The plan of the paper is as follows: §2 introduces the basic formalism 
of both Tree and Particle-Mesh codes. §3 describes the modelling of the 
force in two dimensions. §4 gives the mathematical model for splitting the 
force in two dimensions between the Tree force and the Particle-Mesh force 
components. §5 describes the softening scheme used for the 2d force and 
we analyse errors in force for the 2d TreePM code in §6. We discuss the 
integration of the equations of motion that we use in the 2d TreePM code in 
§7. We also describe a test for self-similar evolution of power law spectra in 
the same section. We present some results of a 2d TreePM simulation run in 
§8. A discussion of the relative merits of the TreePM code and a PM code 
is also given. Computational requirements of our implementation of the 2d 
TreePM code are discussed in §9. 

2. A Review of the Tree and the Particle-Mesh Methods 

2.1 The Tree Method 

We use the same method as the Barnes and Hut (1986) Tree code in our 
implementation of the 2d TreePM code. In a 2d Tree code the simulation 
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area is taken to be a square. If this were to represent the stem of a tree, 
then it will be subdivided at each stage into smaller squares (branches) till 
we reach the particles (leaves). To construct the tree we add particles to the 
simulation area and subdivide any cell that ends up with two particles [6]. 

The force on a particle is computed by adding contribution of other 
particles or of cells. If a cell is too close to the particle, or if it is too big, we 
consider the sub cells of the cell in question instead. The decision is made 
by computing the quantity 9 and comparing it with a threshold Oc- 

e = -<e, (1) 

r 

where d is the size of the cell and r is the distance from the particle to the 
centre of mass of the cell. The error in force increases with Oc- 

The number of terms that contribute to the force of a particle is much 
smaller than the total number of particles for most choices of 6c and this is 
where a tree code gains on a direct force summation method. 

We will use the Barnes and Hut (1986) Tree code, as already mentioned. 
A crucial change to the standard tree walk is that we do not follow nodes 
representing cells that do not have any spatial overlap with the region within 
the threshold radius {vcut^ defined later) for computing the short range force. 

2.2 The Particle-Mesh Method 

A Particle-Mesh (PM) code is the obvious choice for computing long range 
interactions. A PM code adds the construct of a regular grid to the distri- 
bution of particles. The density field represented by particles is interpolated 
onto grid points and the Poisson equation is solved in Fourier space. The 
force is then interpolated back to the positions of particles. Use of a grid 
makes forces inaccurate at the grid scale and smaller scales. In this scheme, 
the mesh and the weight function (Cloud-in-Cell (CIC) in our case) used for 
interpolation between the grid and particle positions are the main sources of 
anisotropy. However, we use the Particle-Mesh method only for computing 
the long range force and errors at small scales do not contribute significantly. 
Also, by de-convolving the interpolating function [4], we reduce errors due 
to anisotropy effects substantially. 

3. The Gravitational Force in Two Dimensions 

When we go from three to two dimensions, we have, in principle, two differ- 
ent ways of modelling the system [2] : (1) We can consider two-dimensional 
perturbations in a three-dimensional expanding Universe. Here we take the 
force between particles to be ^ and assume that all particles (represent- 
ing perturbations) and their velocities are confined to a single plane at the 
initial instant. (2) We can study perturbations that do not depend on one 
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of the three coordinates, i.e., we start with a set of infinitely long straight 
"needles" all pointing along one axis. The force of interaction then falls as 
-. The evolution keeps the "needles" pointed in the same direction, and we 
study the clustering in an orthogonal plane. Particles in the N-body simu- 
lation represent the intersection of these "needles" with this plane. In both 
of these approaches, the Universe is three-dimensional and the background 
is expanding isotropically. Following earlier studies [2], [15], [20], we choose 
the second of the two options. 

More specifically, in order to obtain the force due to perturbations in a 
plane, we solve the Poisson equation for the perturbed part of the gravita- 
tional potential in two dimensions, whereas the unperturbed background is 
still the three dimensional spherically symmetric Friedman Universe. Thus 
the perturbations are described by the mass per unit length, where this 
length is in the direction orthogonal to the two dimensions considered here. 

The gravitational force of a particle situated at the origin in two dimen- 
sions then has the form : 

'Gm' 



f(r) 



r. 



(2) 



Here G is the gravitational coupling constant and m is the mass per unit 
length of the "needle" represented by the particle. 

4. The Mathematical Model for the 2d TreePM Code 



We split the - force into a long range force and a short range force in a 
manner identical to that for the three dimensional TreePM force [3]. We 
compute the long range force in Fourier space and the short range force in 
real space. Following Ewald's method [14], the gravitational potential can 
be split into two parts in Fourier space : 
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where (ff and (/>* are the long range and the short range potentials, respec- 
tively. The splitting is done at the scale r^. 

The expression for the 2d short range force in real space is 



■ exp 
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(6) 
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The above equations describe the mathematical model for force in the 2d 
TreePM code. The long range potential is computed in Fourier space, just as 
in a PM code, but using eqn.4 instead of eqn.3. This potential is then used 
to compute the long range force. The short range force is computed directly 
in real space using eqn.6. This is computed using the Tree approximation. 
The short range force falls rapidly at scales r ':^ Vg and hence we need to 
take this into account only in a small region around each particle. We call 
the scale upto which we add the small scale force as Vcut- 

Evaluation of the short range force can be time consuming. To save time, 
we compute an array containing the magnitude of the short range force at the 
outset. This procedure is identical to that followed in the 3d TreePM code 
[3]. The force between any two objects, particle-cell or particle-particle, is 
then computed by linearly interpolating between the nearby array elements 
multiplied by the unit vector r. It is necessary for the array to sample the 
force at sufficiently closely spaced values of r in order to keep interpolation 
errors in control. 



5. Softening of the Force 

We need to soften the 2d gravitational force at small scales in order to 
ensure coUisionless evolution of the particle distribution in a cosmological 
simulation. We have considered two schemes for softening of the force at 
small scales : 

1) Plummer Softening. The force, in this case, will be given by 

f(r) = - , , ^ ,, r (7) 

where e is the softening length. 

2) Cubic Spline Softening. In this case, we solve the Poisson equation 
in two dimensions with the force due to a point mass replaced by that exerted 
by an extended mass distribution represented by the following : 

p{r) = mW{r,€) (8) 

Here W{r, e) is the normalised spline kernel used in the SPH formalism [19] 
with e the smoothing length. 

W{r, e) has the following form in two dimensions. 

/ 40 \ f l-6(7)' + 6(i)', 0>i<0.5 
^^"''^=7^ 2(l-f)^ 0.5<i<1.0 (9) 

^ ^ [ f >1.0 

Solving the Poisson equation and using relevant boundary conditions for 
the potential and its first derivative (i.e. the force) to obtain the constants 
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of integration, we get the cubic spline softened potential and then obtain 
the force f as a gradient of the potential. 

> f < 0.5; 
f(r) = <' _(80_160r 12(>£_32ri_^|^r 0.5 < ^ < 1.0; (10) 

f >1.0 

One can easily see that the spline softening has the advantage that the 
force becomes exactly Newtonian for r > e, while the Plummer force con- 
verges relatively slowly to the Newtonian form. Dehnen [11] has argued that 
compact softening kernels are superior and we use the spline softened kernel 
in our implemention of the 2d TreePM code. 



6. Error Estimation 

It is important to estimate the errors in numerical evaluation of force in a 
realistic situation, even though we do not expect errors to add up coher- 
ently. For a comprehensive study in errors in force introduced by various 
components of the TreePM code we would like to refer to [4]. Though the 
above-mentioned work is in the context of the 3d TreePM code, the key 
features of the analysis actually carry over to the 2d TreePM code as well. 
We test errors for two distributions of particles : a homogeneous distri- 
bution and a clumpy distribution. For the homogeneous distribution, we use 
randomly distributed particles in a box. We use 1024^ particles distributed 
on a 1024^ grid. We compute the force using a reference setup (rg = 4, 
6c = 0.01, rcut = 6rs) and the setup we wish to test (r^ = 1, 9c = 0.5, 
i^cut = 5rs). We compute the fractional error in force acting on each parti- 
cle. This is defined as 

_ |f — fre/l 



Ire/ 1 



(11) 



Fig.l shows the cumulative distribution of fractional errors. The curves 
show the fraction of particles with error greater than e. The thick line 
shows this for the homogeneous distribution. Error e for 99% of particles 
is less than 4%. Results for the clumpy distribution of particles are shown 
by the dashed line. For this case, we used the output of a 2d power law 
(n = 1) simulation run in an Einstein deSitter background Universe with 
the TreePM code. Errors in this case are somewhat smaller as compared 
to the homogeneous distribution for much the same reason as that for a 3d 
Tree code [17] or a 3d TreePM code [3]. Error e for 99% of particles is less 
than 2% for the clumpy distribution. 
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Figure 1. This figure shows the distribution of errors. The variation of the fraction 
of particles with error greater than a threshold, as a function of the threshold error, 
is plotted. Thick line marks the error for a homogeneous distribution of particles 
and the dashed line shows the same for a clumpy distribution. These errors were 
measured with respect to a reference force, determined with a very conservative 
value of Ts and 9^- This panel shows that 99% of the particles have fractional error 
in force that is less than 4% for the homogeneous distribution and less than 2% for 
the clumpy distribution. 



7. Integrating the Equation of Motion 

Our discussion so far has dealt only with the evaluation of force. This is 
the main focus of this paper as the key difference between the TreePM and 
other methods is in the scheme used for evaluation of force. However for the 
sake of completeness, we give here details of integration of the equations of 
motion used in the code. We use an Einstein deSitter background cosmology 
for all our 2d simulations. The equations of motion are then given by the 
simple form 



X + 2-X 
a 






8 Suryadeep Ray 

V^ = A7TGa'^{p-p) (12) 

Here x is the comoving coordinate, a is the scale factor, (p is the gravitational 
potential of perturbations, p is the total density and p is the average density 
of the Universe. Dot represents differentiation with respect to time. We can 
recast these equations in the following form. 

x" + — x' = -— W 
2a 2a 

VV = (5 = --! (13) 

P 

Here prime denotes differentiation with respect to the scale factor, S is the 
density contrast and ip is the appropriately scaled gravitational potential of 
perturbations. 

The equations of motion are identical to that in three dimensions apart 
from the fact that all the vectors under consideration are two-dimensional 
vectors. The functional form of the gravitational force given by V0 is, of 
course, different. 

One can see that the equations of motion contain a velocity dependent 
term and hence we cannot use the usual leap-frog method. We recast the 
leap-frog method so that velocities and positions are defined at the same 
instant [4]. We solve the equation for velocity iteratively. Time step is 
chosen to be a small fraction of the smallest dynamical time in the system 
at any given stage. The fraction to be chosen is fixed by checking for scale 
invariance in evolution of power law spectra : a simulation is repeated with 
different choices of timestep until we find the largest timestep for which we 
can reach the highly non-linear regime and retain scale invariance as well. 
We then use a timestep that is half of this largest time step. 

Fig. 2 shows ^ as a function of r/rni{t) for several epochs obatined from 
a 2d TreePM simulation of a power law model with index n = —0.4. r„;(t) 
is the scale which is going non-linear at time t and it varies in proportion 
with o^' '""''^^ in the Einstein deSitter model. In the figure, we have only 
plotted ^ at scales more than two times larger than the artificial softening 
length used in the simulation. We can see that scale invariance holds for the 
spectrum over a wide range which means that we can probe the non-linear 
regime in gravitational clustering with a high degree of accuracy using the 
2d TreePM code. 

8. The 2d TreePM Code vs. the 2d Particle-Mesh Code 

In this section, we present a brief comparison of the 2d TreePM and Particle- 
Mesh methods with the aim of highlighting the efficacy of our method in 2d 
cosmological simulations. 

We ran a 2d simulation of a power law model with index n = —0.4 with 
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Figure 2. This figure shows ^ as a function of r/rni{t) for several epochs. Here 
rni{t) is the scale which is going non-linear at time t and it varies in proportion 
with a^/(^+^> in the Einstein deSitter model. The index of the power spectrum is 
n = —0.4. We have only plotted ^ at scales more than two times larger than the 
artificial softening length used in the simulation. 



1024^ particles on a 1024^ grid in an Einstein deSitter background Universe 
with a PM code as well as with the TreePM code discussed here. For the 
TreePM run we used rg = 1, 9c = 0.5 and softening parameter e = 0.2. 

The top panel of Fig. 3 shows identical boxes from two independent sim- 
ulations with the same initial conditions. The top left panel shows a sim- 
ulation with the TreePM code and the top right panel shows the same for 
a PM code. The large scale structures appear to be the same in the two. 
However, one can see that there are significant differences at small scales 
when one plots the two-point correlation function for the two cases, ^(r) 
is plotted as a function of scale r in the bottom left panel of fig. 3. The 
thick line shows the correlation function for the TreePM simulation and the 
dashed line shows the same for the PM simulation. We have only plotted ^ 
at scales more than two times larger than the artificial softening length used 
in the TreePM simulation. The correlation function in the TreePM simula- 
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Figure 3. The top panel of this figure shows a box from a simulation of a power 
law model with index n = —0.4 for two cases. The top left panel shows the box from 
a TreePM simulation. For comparison, we have included the same box from a PM 
simulation of the same initial conditions in the top right panel. For convenience, we 
have randomly plotted only one out of every sixteen particles in the original simula- 
tions in the figure. The large scale structures appear to be the same in the two. We 
will look more closely at the circled haloes in both panels in the following discus- 
sion. The bottom left panel of this figure shows the averaged correlation function 
S,(r) as a function of scale in grid units. The thick line shows this quantity for the 
TreePM simulation and the dashed line shows the same for the PM simulation. We 
have only plotted ^ at scales more than two times larger than the artificial softening 
length used in the TreePM simulation. The point marked by an arrow in the figure 
represents the larger softening scale for the PM code. The correlation functions 
match at large scales but the PM simulation underestimates the clustering at small 
scales. The bottom right panel is a plot which shows average density p within a 
sphere of radius r as a function of r (in grid units) for the two halos circled in the 
top panel of this figure. Again, the thick line shows average density for the TreePM 
simulation and the dashed line shows the same for the PM simulation. Here also we 
have plotted average density at scales more than two times larger than the artificial 
softening length used in the TreePM simulation. The density profiles match at large 
scales as expected, but one can see that the TreePM simulation gives rise to haloes 
with higher central densities. 
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tion matches with that from the PM simulation at large scales, but at scales 
of the order of unity (in grid units) and below, the TreePM simulation has 
a higher correlation function. This is to be expected because of the superior 
force resolution of the TreePM method as opposed to the PM force, where 
the force is softened naturally at the grid scale. The scale of softening for 
the PM code is marked by an arrow in the figure. 

We also study the density profiles of the two circled haloes in the top 
panel of fig. 3. These particular haloes have been chosen as representatives 
for the analysis because they are reasonably large, dense and spherically 
symmetric. The bottom right panel of fig. 3 shows average density p within 
a sphere of radius r from the halo centre plotted as a function of r for the 
two haloes. The full line shows the density profile for the halo from the 
TreePM simulation and the dashed line the same from the PM simulation. 
Here we have only plotted average density at scales more than four times 
larger than the artificial softening length used in the TreePM simulation. 
We can see that, though not visibly obvious from fig. 3, the halo from the 
TreePM simulation is clearly far denser in the central region as compared to 
the halo from the Particle-Mesh simulation. The density profiles converge 
at some distance from the halo centres as expected. 

9. Computational Requirements 

In this section, we describe the computational resources required for the 
present implementation of the 2d TreePM code. Given that we have com- 
bined the Tree and the PM codes, the memory requirement is obviously 
greater than that for either one code. We need four arrays for the PM part 
i.e. for the potential and the force. The rest is exactly the same as a stan- 
dard Barnes and Hut 2d Tree code. With efficient memory management, 
we need less than 75MB of RAM for a simulation with 1024^ particles on a 
1024^ grid. The number mentioned is for floating point variables. If we use 
double precision variables, our requirement will go up by a factor of two. 
The time taken (per time step per particle) by the 2d TreePM code 

{rs = 1, ^c = 0.5, rent = 4.5r„ Npartide = 10242, Ngrid = 10242) is of the 

order of 240 microseconds. This number was generated using a 2.4GHz Xeon 
personal computer where the code was compiled with the Intel F90 compiler. 

10. Discussion 

In this paper, we have described the two-dimensional TreePM method in 
detail. Our method offers greater dynamic range and superior resolution as 
compared to a 2d Particle-Mesh method and can therefore probe the non- 
linear regime in two-dimensional cosmological simulations more effectively. 
We believe that a 2d TreePM code will allow us to explore a higher dynamic 
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range in densities (and ^) for studying scaling relations in two-dimensions 
as compared to earlier work done using Particle-Mesh codes [2]. Work is in 
progress in this direction and will be reported elsewhere. 

The 2d TreePM code is also amenable to parallelisation along the lines 
of the 3d TreePM code [5, 23] and is likely to scale well. 
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